knitr::opts_chunk$set(fig.width=5, fig.height=5) 

rm(list=ls())


options(scipen=999)

pkgs <- c("dplyr", "ggplot2", "stargazer", "car", "tidyverse", "extrafont", "haven", "coefplot", "plotrix", "sensemakr")
# A function to load all above packages. Install if they have not been installed.
usePackage <- function(p){
  for (pkg in p){
    if (!is.element(pkg, installed.packages()[,1]))
      install.packages(pkg, dep = TRUE, repos = "https://cloud.r-project.org/")
    require(pkg, character.only = TRUE)
  }
}
usePackage(pkgs)


d1 <- read.csv("Study 2 Full Data.csv")


d1$Duplicate <- ifelse(d1$Q_RelevantIDDuplicate == "TRUE", 1, 0)
d1$Bot <- ifelse(d1$Q_RecaptchaScore == "NA", 1, 0 )

d0 <- subset(d1, is.na(Duplicate) & d1$Bot == 0)



d0$Republican <- 0
d0$Republican[d0$PartyID == "Republican" | d0$PartyID2 == "Closer to Republican Party"] <- 1
d0$Democrat <- 0
d0$Democrat[d0$PartyID == "Democrat" | d0$PartyID2 == "Closer to Democratic Party"] <- 1

d0$Independent <- as.numeric(ifelse(d0$Democrat == 1 | d0$Republican == 1, 0, 1))

d0$Legislator <- ifelse(d0$Job == "State legislator", 1, 0)


d <- subset(d0, d0$Independent == 0 & d0$Legislator == 1)

write.csv(d, "Study 2 Clean Data.csv")

d$FTDemVoters <- d$FT_1
d$FTRepVoters <- d$FT_2
d$FTDemPols <- d$FT_5
d$FTRepPols <- d$FT_6

d$FTOutPartyPeople <- ifelse(d$Republican == 1, d$FTDemVoters, d$FTRepVoters)
d$FTInPartyPeople <- ifelse(d$Republican == 1, d$FTRepVoters, d$FTDemVoters)


d$FTOutPartyPols <- ifelse(d$Republican == 1, d$FTDemPols, d$FTRepPols)
d$FTInPartyPols <- ifelse(d$Republican == 1, d$FTRepPols, d$FTDemPols)

d$AffPolPeople <- d$FTInPartyPeople - d$FTOutPartyPeople


d$AffPolPols <- d$FTInPartyPols - d$FTOutPartyPols

d$AffPol <- (d$AffPolPols + d$AffPolPeople)/2


## AFF POL GRAPHS:

library(dplyr)
library(tidyr)
library(ggplot2)

d <- d %>%
  mutate(
    FTOutPartyPeople = ifelse(Republican == 1, FTDemVoters, FTRepVoters),
    FTInPartyPeople = ifelse(Republican == 1, FTRepVoters, FTDemVoters),
    FTOutPartyPols = ifelse(Republican == 1, FTDemPols, FTRepPols),
    FTInPartyPols = ifelse(Republican == 1, FTRepPols, FTDemPols)
  )

# Reshape data for plotting
d_long <- d %>%
  select(FTOutPartyPeople, FTInPartyPeople, FTOutPartyPols, FTInPartyPols) %>%
  pivot_longer(cols = everything(), names_to = "Category", values_to = "FeelingThermometer") %>%
  mutate(
    PartyType = case_when(
      str_detect(Category, "OutParty") ~ "Out-party",
      str_detect(Category, "InParty") ~ "In-party"
    ),
    Group = case_when(
      str_detect(Category, "People") ~ "People",
      str_detect(Category, "Pols") ~ "Politicians"
    )
  ) %>%
  filter(!is.na(FeelingThermometer) & is.finite(FeelingThermometer))


# Create the density plot with HEX colors
plot <- ggplot(d_long, aes(x = FeelingThermometer, fill = PartyType, color = PartyType)) +
  geom_density(alpha = 0.6) +
  facet_wrap(~ Group, scales = "fixed") +  # Ensuring the y-axis is the same for both facets
  scale_fill_manual(name = "Party Type", values = c("Out-party" = "#FF7F50", "In-party" = "#008080")) +
  scale_color_manual(name = "Party Type", values = c("Out-party" = "#FF7F50", "In-party" = "#008080")) +
  theme_minimal() +
  theme(
    text = element_text(family = "LM Roman 10"),
    plot.title = element_text(hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill = NA, size = .5),
    legend.position = "bottom",
    strip.background = element_rect(fill = "gray90", color = "black", size = 0.5),
    strip.text = element_text(face = "bold")
  ) +
  labs(
    x = "Feeling Thermometer",
    y = "Proportion of Respondents",
    title = "In-party vs. Out-party Feelings: People vs. Politicians",
    fill = "Party Type"
  )

# Display the plot
print(plot)

# Save the plot
ggsave(filename = "In_Out_Party_Feelings_Density_FixedY_NewColors.jpeg", plot = plot, device = "jpeg", height = 6, width = 8, units = "in", dpi = 500)



d$RedistTreatment <- as.factor(ifelse(d$DistrictCondition == "DistrictInParty", 1, 0))
d$JudgesTreatment <- as.factor(ifelse(d$JudgesCondition == "JudgesInParty", 1, 0))
d$JudgesTreatment2 <- ifelse(d$JudgesTreatment == 1, "Yes", "No")


d$JudgesResponse0 <- ifelse(d$JudgesTreatment == 1, d$JudgesInParty, d$JudgesControl)
d$JudgesResponse <- case_when(
  d$JudgesResponse0 == "Strongly support" ~ 3,
  d$JudgesResponse0 == "Support" ~ 2,
  d$JudgesResponse0 == "Slightly support" ~ 1,
  d$JudgesResponse0 == "Neither support nor oppose" ~ 0,
  d$JudgesResponse0 == "Slightly oppose" ~ -1,
  d$JudgesResponse0 == "Oppose" ~ -2,
  d$JudgesResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$JudgesResponse <- (d$JudgesResponse + 3) / 6
d$JudgesResponse <- as.numeric(d$JudgesResponse)

summary(lm(d$JudgesResponse ~ d$JudgesTreatment))

m1 <- lm(d$JudgesResponse ~ d$JudgesTreatment)

## SECOND ORDER BELIEFS AND PLOTS

d$RedistResponse0 <- ifelse(d$RedistTreatment == 1, d$DistrictInParty, d$DistrictControl)
d$RedistResponse <- case_when(
  d$RedistResponse0 == "Strongly support" ~ 3,
  d$RedistResponse0 == "Support" ~ 2,
  d$RedistResponse0 == "Slightly support" ~ 1,
  d$RedistResponse0 == "Neither support nor oppose" ~ 0,
  d$RedistResponse0 == "Slightly oppose" ~ -1,
  d$RedistResponse0 == "Oppose" ~ -2,
  d$RedistResponse0 == "Strongly oppose" ~ -3,
  TRUE ~ NA_real_
)
d$RedistResponse <- (d$RedistResponse + 3) / 6
d$RedistResponse <- as.numeric(d$RedistResponse)

d$RedistTreatment2 <- ifelse(d$RedistTreatment == 1, "Yes", "No")

summary(lm(d$RedistResponse ~ d$RedistTreatment  ))

m2 <- lm(d$RedistResponse ~ d$RedistTreatment)

## Beliefs about others - DISTRICTS

d$SORedistDems <- case_when(
  d$DistrictDems == "All" ~ 3,
  d$DistrictDems == "Almost all" ~ 2,
  d$DistrictDems == "Most" ~ 1,
  d$DistrictDems == "About half" ~ 0,
  d$DistrictDems == "Less than half" ~ -1,
  d$DistrictDems == "Almost none" ~ -2,
  d$DistrictDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SORedistDems <- (d$SORedistDems + 3) / 6
d$SORedistDems <- as.numeric(d$SORedistDems)

#summary(lm(d$SORedistDems ~ d$Democrat  ))

soDm1 <- lm(d$SORedistDems ~ d$Democrat  )


d$SORedistReps <- case_when(
  d$DistrictReps == "All" ~ 3,
  d$DistrictReps == "Almost all" ~ 2,
  d$DistrictReps == "Most" ~ 1,
  d$DistrictReps == "About half" ~ 0,
  d$DistrictReps == "Less than half" ~ -1,
  d$DistrictReps == "Almost none" ~ -2,
  d$DistrictReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SORedistReps <- (d$SORedistReps + 3) / 6
d$SORedistReps <- as.numeric(d$SORedistReps)

soRm1 <- lm(d$SORedistReps ~  d$Republican )

#summary(lm(d$SORedistReps ~  d$Republican ))


d$SORedistInParty <- ifelse(d$Democrat == 1, d$SORedistDems, d$SORedistReps)
d$SORedistOutParty <- ifelse(d$Democrat == 1, d$SORedistReps, d$SORedistDems)


d$RedistTreatmentResponse <- as.numeric(ifelse(d$RedistTreatment == 1, d$RedistResponse, "NA"))
d$RedistControlResponse <- as.numeric(ifelse(d$RedistTreatment == 0, d$RedistResponse, "NA"))


#summary(d$RedistTreatmentResponse)
#summary(d$SORedistInParty)
#summary(d$SORedistOutParty)



mean_data <- data.frame(
  variable = factor(c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support"), levels = c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support")[order(c(mean(d$RedistControlResponse, na.rm = TRUE), mean(d$RedistTreatmentResponse, na.rm = TRUE), mean(d$SORedistOutParty, na.rm = TRUE), mean(d$SORedistInParty, na.rm = TRUE)))]),
  mean = c(mean(d$RedistControlResponse, na.rm = TRUE),
           mean(d$RedistTreatmentResponse, na.rm = TRUE), 
           mean(d$SORedistOutParty, na.rm = TRUE), 
           mean(d$SORedistInParty, na.rm = TRUE)),
  se = c(sd(d$RedistControlResponse, na.rm = TRUE) / sqrt(sum(!is.na(d$RedistControlResponse))),
         sd(d$RedistTreatmentResponse, na.rm = TRUE) / sqrt(sum(!is.na(d$RedistTreatmentResponse))),
         sd(d$SORedistOutParty, na.rm = TRUE) / sqrt(sum(!is.na(d$SORedistOutParty))),
         sd(d$SORedistInParty, na.rm = TRUE) / sqrt(sum(!is.na(d$SORedistInParty))))
)

soPA1 <- ggplot(mean_data, aes(x = variable, y = mean, fill = variable)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + geom_hline(yintercept = 0.5, linetype = "dashed", color = "black", position = "identity") +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "Scrap bipartisan committee and allow governor to appoint\n all members of redistricting committee",
       x = "", y = "Mean Response") +
  scale_x_discrete(labels = c("Control - Support\n proposal", "Treatment - Support\n in-party proposal", "Belief about\n in-party legislator\n support for in-party\n proposal", "Belief about\n out-party legislator\n support for out-party\n proposal")) +  theme( axis.line.y = element_blank(), axis.ticks.y=element_blank())+   theme(text         = element_text(size=10, family="LM Roman 10")) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +theme(plot.title = element_text(hjust = 0.5, family = "LM Roman 10"))  +theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position = "none") +scale_y_continuous( labels = scales::percent_format(accuracy = 1),  limits =  c(0, 1))

#soPA1

#ggsave(filename="soP1.jpeg", plot=soPA1, device="jpeg", height=5, width=5, units="in", dpi=500)



## Beliefs about others - JUDGES

d$SOJudgesDems <- case_when(
  d$JudgesDems == "All" ~ 3,
  d$JudgesDems == "Almost all" ~ 2,
  d$JudgesDems == "Most" ~ 1,
  d$JudgesDems == "About half" ~ 0,
  d$JudgesDems == "Less than half" ~ -1,
  d$JudgesDems == "Almost none" ~ -2,
  d$JudgesDems == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOJudgesDems <- (d$SOJudgesDems + 3) / 6
d$SOJudgesDems <- as.numeric(d$SOJudgesDems)

#summary(lm(d$SOJudgesDems ~ d$Democrat  ))
soDm2 <- lm(d$SOJudgesDems ~ d$Democrat  )

d$SOJudgesReps <- case_when(
  d$JudgesReps == "All" ~ 3,
  d$JudgesReps == "Almost all" ~ 2,
  d$JudgesReps == "Most" ~ 1,
  d$JudgesReps == "About half" ~ 0,
  d$JudgesReps == "Less than half" ~ -1,
  d$JudgesReps == "Almost none" ~ -2,
  d$JudgesReps == "None" ~ -3,
  TRUE ~ NA_real_
)
d$SOJudgesReps <- (d$SOJudgesReps + 3) / 6
d$SOJudgesReps <- as.numeric(d$SOJudgesReps)

#summary(lm(d$SOJudgesReps ~   d$Republican ))
soRm2 <- lm(d$SOJudgesReps ~   d$Republican )

d$SOJudgesInParty <- ifelse(d$Democrat == 1, d$SOJudgesDems, d$SOJudgesReps)
d$SOJudgesOutParty <- ifelse(d$Democrat == 1, d$SOJudgesReps, d$SOJudgesDems)


d$JudgesTreatmentResponse <- as.numeric(ifelse(d$JudgesTreatment == 1, d$JudgesResponse, "NA"))
d$JudgesControlResponse <- as.numeric(ifelse(d$JudgesTreatment == 1, d$JudgesResponse, "NA"))

#summary(d$JudgesTreatmentResponse)
#summary(d$SOJudgesInParty)
#summary(d$SOJudgesOutParty)

mean_data2 <- data.frame(
  variable = factor(c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support"), levels = c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support")[order(c(mean(d$JudgesControlResponse, na.rm = TRUE), mean(d$JudgesTreatmentResponse, na.rm = TRUE), mean(d$SOJudgesOutParty, na.rm = TRUE), mean(d$SOJudgesInParty, na.rm = TRUE)))]),
  mean = c(mean(d$RedistControlResponse, na.rm = TRUE),
           mean(d$JudgesTreatmentResponse, na.rm = TRUE), 
           mean(d$SOJudgesOutParty, na.rm = TRUE), 
           mean(d$SOJudgesInParty, na.rm = TRUE)),
  se = c(sd(d$JudgesControlResponse, na.rm = TRUE) / sqrt(sum(!is.na(d$JudgesControlResponse))),
         sd(d$JudgesTreatmentResponse, na.rm = TRUE) / sqrt(sum(!is.na(d$JudgesTreatmentResponse))),
         sd(d$SOJudgesOutParty, na.rm = TRUE) / sqrt(sum(!is.na(d$SOJudgesOutParty))),
         sd(d$SOJudgesInParty, na.rm = TRUE) / sqrt(sum(!is.na(d$SOJudgesInParty))))
)

soPA2 <- ggplot(mean_data2, aes(x = variable, y = mean, fill = variable)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + geom_hline(yintercept = 0.5, linetype = "dashed", color = "black", position = "identity") +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "Appoint judges in groups to reduce minority party\n ability to block individual judges",
       x = "", y = "Mean Response") + scale_x_discrete(labels = c("Control - Support\n proposal", "Treatment - Support\n in-party proposal", "Belief about\n in-party legislator\n support for in-party\n proposal", "Belief about\n out-party legislator\n support for out-party\n proposal")) +  theme( axis.line.y = element_blank(), axis.ticks.y=element_blank())+   theme(text         = element_text(size=10, family="LM Roman 10")) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +theme(plot.title = element_text(hjust = 0.5, family = "LM Roman 10"))  +theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position = "none") +scale_y_continuous( labels = scales::percent_format(accuracy = 1),  limits =  c(0, 1))

#soPA2

#ggsave(filename="soP2.jpeg", plot=soPA2, device="jpeg", height=5, width=5, units="in", dpi=500)


## MAKE THE DATA LONG 


##Response_Type variables 
d$JResponse <- d$JudgesResponse
d$RResponse <- d$RedistResponse

## Treatment Type variables
d$TreatJ <- d$JudgesTreatment
d$TreatR <- d$RedistTreatment



y <- dplyr::mutate(d, ID = row_number())


## DVs to common stem
y$DV_1 <- y$RResponse
y$DV_2 <- y$JResponse
y$DV_3 <- y$SOJudgesInParty
y$DV_4 <- y$SOJudgesOutParty
y$DV_5 <- y$SORedistInParty
y$DV_6 <- y$SORedistOutParty

## Reshape to long
longS <- reshape(y, direction = "long", 
                 varying = list(c("DV_1","DV_2","DV_3","DV_4","DV_5","DV_6")))
longS$DV <- longS$time
longS$response <- longS$DV_1

## Code treatment variable
longS$treat <- ifelse(longS$DV == 1, as.numeric(longS$TreatR), 
                      ifelse(longS$DV == 2, as.numeric(longS$TreatJ), NA)
                      )
longS$treat <- longS$treat - 1
table(longS$treat)


## Code in vs out party
longS$inorout <- ifelse(longS$DV %in% c(3,5), 1, 
                        ifelse(longS$DV %in% c(4,6), 0, NA))

## Code judge or redistrict
longS$JorR <- ifelse(longS$DV %in% c(3,4), 1, 
                        ifelse(longS$DV %in% c(5,6), 0, NA))

## Model for experiment

m_1 <- lm(response ~ treat + I(DV - 1), data = subset(longS, DV == 1 | DV == 2))
summary(m_1)

# clustered SEs
library(sandwich)
ses_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% 1:2])))
t_1 <- abs(coef(m_1) / ses_1)
round(cbind(coef(m_1), 
      ses_1,
      t_1,
      (1 - pt(t_1, m_1$df.residual))*2
       ),
      3)

coefs_clustered_m1 <- cbind(coef(m_1), ses_1, t_1, (1 - pt(t_1, m_1$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered_m1) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered_m1, 
          title = "Model 1", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))



## Model for SO beliefs

m_2 <- lm(response ~ inorout*AffPol + JorR, data = subset(longS, DV %in% 3:6))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)

coefs_clustered <- cbind(coef(m_2), ses_2, t_2, (1 - pt(t_2, m_2$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered, 
          title = "SO Beliefs Model", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))


control_mean_1 <- coef(m_1)[1]
control_se_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% 1:2])))[1]

treat_mean_1 <- coef(m_1)[1] + coef(m_1)[2]
treat_se_1 <- sqrt(control_se_1^2 + (vcovCL(m_1, cluster = longS$id[longS$DV %in% 1:2]))[1, 2]^2)

# Calculate mean and standard errors for intercept and "inorout" variable in m_2 model
control_mean_2 <- coef(m_2)[1]
control_se_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))[1]

inorout_mean_2 <- coef(m_2)[1] + coef(m_2)[2]
inorout_se_2 <- sqrt(control_se_2^2 + (vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6]))[1, 2]^2)

# Create table
table_data <- data.frame(
  variable = c("Control", "Support the proposal", "Beliefs about in-party support", "Beliefs about out-party support"),
  mean = c(control_mean_1, treat_mean_1, inorout_mean_2, control_mean_2),
  se = c(control_se_1, treat_se_1, inorout_se_2, control_se_2)
)

mean_data3 <- as.data.frame(table_data)
mean_data3$variable <- factor(mean_data3$variable, levels = c("Control","Support the proposal", "Beliefs about in-party support",  "Beliefs about out-party support"))




soplot <- ggplot(mean_data3, aes(x = variable, y = mean, fill = variable)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) + geom_hline(yintercept = 0.5, linetype = "dashed", color = "black", position = "identity") +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "State legislator support for policy proposals\n that would reduce the power of the party in the minority",
       x = "", y = "Mean Response") + scale_x_discrete(labels = c("Control:\n Non-partisan Support", "Treatment:\n Partisan Support", "Perceived\n In-Party Support", "Perceived\n Out-Party Support")) +  theme( axis.line.y = element_blank(), axis.ticks.y=element_blank())+   theme(text         = element_text(size=10, family="LM Roman 10")) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +theme(plot.title = element_text(hjust = 0.5, family = "LM Roman 10"))  +theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+ theme(legend.position = "none") +scale_y_continuous( labels = scales::percent_format(accuracy = 1),  limits =  c(0, 1))

soplot


ggsave(filename="soplot.jpeg", plot=soplot, device="jpeg", height=5.5, width=5.5, units="in", dpi=500)

mean_data4 <- head(mean_data3, nrow(mean_data3)-2)

soplotExp <- ggplot(mean_data4, aes(x = variable, y = mean, fill = variable)) +
  geom_point(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_errorbar(aes(ymin = mean - 1.96 * se, ymax = mean + 1.96 * se),
                width = 0.2, position = position_dodge(width = 0.9)) +
  labs(title = "State legislator support for policy proposals\n that would reduce the power of the party in the minority",
       x = "",
       y = "Mean level of support (50% = neither support nor oppose)") +
  scale_x_discrete(labels = c("Control -\n Support for proposal\n (No partisan information)", 
                              "Treatment -\n Support for proposal\n (Benefits in-party)")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 1)) +
  theme_bw() +
  theme(axis.line.y = element_blank(),
        axis.ticks.y = element_blank(),
        text = element_text(size=10, family="LM Roman 10"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"),
        legend.position = "none")

soplotExp

ggsave(filename="soplotExp.jpeg", plot=soplotExp, device="jpeg", height=5, width=5, units="in", dpi=500)

## Partisan heterogeneity for Table 7


m_1 <- lm(response ~ treat*Republican  + I(DV - 1), data = subset(longS, DV == 1 | DV == 2))
summary(m_1)

# clustered SEs
library(sandwich)
ses_1 <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% 1:2])))
t_1 <- abs(coef(m_1) / ses_1)
round(cbind(coef(m_1), 
            ses_1,
            t_1,
            (1 - pt(t_1, m_1$df.residual))*2
),
3)

coefs_clustered_m1_updated <- cbind(coef(m_1), ses_1, t_1, (1 - pt(t_1, m_1$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered_m1_updated) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered_m1_updated, 
          title = "Model 1 (Updated)", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))


library(ggplot2)
library(dplyr)
library(broom)

# Calculate coefficients and their clustered standard errors
tidy_model_1 <- tidy(m_1, conf.int = TRUE)

# Calculate the clustered standard errors
cluster_se <- sqrt(diag(vcovCL(m_1, cluster = longS$id[longS$DV %in% c(1:2)])))

# Calculate the t values
t_values <- abs(tidy_model_1$estimate / cluster_se)

# Calculate the p values
p_values <- 2 * (1 - pt(t_values, df = m_1$df.residual))

tidy_model_1 <- tidy_model_1 %>%
  mutate(`Clustered SE` = cluster_se,
         `t value` = t_values,
         `p value` = p_values)

# Rename variables for more intuitive plotting and filter out "(Intercept)" and "I(DV - 1)" coefficients
tidy_model_1 <- tidy_model_1 %>%
  rename(Coefficient = estimate, term = term) %>%
  filter(!term %in% c("(Intercept)", "I(DV - 1)")) %>%
  mutate(term = factor(term, levels = c("treat", "Republican", "AffPolS"),
                       labels = c("Treatment\n effect", "Republican\n identification", 
                                  "Affective\n polarization"))) %>%
  mutate(conf.low.90 = Coefficient - qnorm(0.95)*`Clustered SE`,
         conf.high.90 = Coefficient + qnorm(0.95)*`Clustered SE`,
         conf.low.95 = Coefficient - qnorm(0.975)*`Clustered SE`,
         conf.high.95 = Coefficient + qnorm(0.975)*`Clustered SE`)
# Generate the coefficient plot with 90% and 95% confidence intervals
cplot1 <- ggplot(tidy_model_1, aes(y = term, x = Coefficient)) +
  geom_vline(xintercept = 0, color = "black") +
  geom_point(color = "black") +
  geom_segment(aes(x = conf.low.90, xend = conf.high.90, y = term, yend = term), size = 1, color = "black") +
  geom_segment(aes(x = conf.low.95, xend = conf.high.95, y = term, yend = term), size = 0.5, color = "black") +
  labs(title = "Percent increase in\n support for policies\n subverting minority power", 
       x = "", 
       y = "") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 2), limits = c(-.1, .3)) +
  theme_bw() +
  theme(text = element_text(size = 10, family = "LM Roman 10"),
        plot.title = element_text(hjust = 0.5, family = "LM Roman 10"),
        axis.line.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"),
        legend.position = "none")

ggsave(filename="coefplot.jpeg", plot=cplot1, device="jpeg", height=5, width=5, units="in", dpi=500)

longS$outorin <-  ifelse(longS$inorout == 0, 1, 0)

m_2 <- lm(response ~ outorin + Republican + AffPolS + JorR, data = subset(longS, DV %in% 3:6))
summary(m_2)

# clustered SEs
ses_2 <- sqrt(diag(vcovCL(m_2, cluster = longS$id[longS$DV %in% 3:6])))
t_2 <- abs(coef(m_2) / ses_2)
round(cbind(coef(m_2), 
            ses_2,
            t_2,
            (1 - pt(t_2, m_2$df.residual))*2
),
3)

coefs_clustered_m2 <- cbind(coef(m_2), ses_2, t_2, (1 - pt(t_2, m_2$df.residual)) * 2)

# Name the columns
colnames(coefs_clustered_m2) <- c("Estimate", "Clustered SE", "t-value", "p-value")

# Convert the data frame to a LaTeX table using stargazer
stargazer(coefs_clustered_m2, 
          title = "Model 2", 
          summary = FALSE, 
          header = FALSE, 
          type = "latex", 
          digits = 3,
          column.labels = c("Estimate", "Clustered SE", "t-value", "p-value"))


#cplot2

#ggsave(filename="coefplot2.jpeg", plot=cplot2, device="jpeg", height=6, width=6, units="in", dpi=500)

## Test for Partisan :


